Multi-OMICs Factor Analysis (MOFA) and Unsupervised Data Learning

Multi-OMICs Factor Analysis (MOFA)

MOFA is a typical hypothesis-free data exploration framework, https://www.embopress.org/doi/10.15252/msb.20178124. It allows data integration by extracting common axes of variation from the multiple OMICs layers. Given several data matrices with measurements of multiple ‘omics data types on the same or on overlapping sets of samples, MOFA infers a low-dimensional data representation in terms of (hidden) factors which can be then visualized and interpreted.

Importantly, MOFA disentangles whether the underlying axes of heterogeneity are unique to a single data modality or are manifested in multiple modalities, thereby identifying links between the different ‘omics. Once trained, the model output can be used for a range of downstream analyses, including

• Visualisation of samples in factor space
• Annotation of factors using (gene set) enrichment analysis
• Imputation of missing values
• Support of sparse data
• Support of OMICs with non-Gaussian distribution including binary and count data

How does a typical MOFA output look like? In the MOFA paper they applied this method to Chronic Lymphocytic Leukaemia (CLL) on 200 human patients that combined 1) drug response, 2) somatic mutations (targeted sequencing), 3) RNAseq data, and 4) Methylation array data.

Notably, nearly 40% of the 200 samples were profiled with some but not all OMICs types; such a missing value scenario is not uncommon in large cohort studies, and MOFA is designed to cope with it.

One can inspect loadings of MOFA hidden factors and observe known bio-markers fo CLL. Specifically, factor 1 was aligned with the somatic mutation status of the immunoglobulin heavy-chain variable region gene (IGHV), while Factor 2 aligned with trisomy of chromosome 12.

Loadings from each factor can be used for pathway enrichment analysis. Interestingly, markers from diffeerent OMICs co-varying in a certain MOFA factor can be shown to be coupled to the same biological pathway. For CLL data, factor 5 demonstrates co-variation between drug response and mRNA levels. This factor tagged a set of genes that were highly enriched for oxidative stress and senescence pathways. Consistent with the annotation based on the mRNA view, it is observed that the drugs with the strongest weights on Factor 5 were associated with response to oxidative stress, such as target reactive oxygen species, DNA damage response and apoptosis.

Next MOFA can perform efficient imputation of missing values including imputation of "missing patients", i.e. when a patient is present in one OMICs layer but absent in another, this absent layer can be reconstructed for the patient, i.e. gray bars in the panel a) of the figure above can be filled.

Finally, latent factors inferred by MOFA are predictive of clinical outcomes. The authors of MOFa paper assessed the prediction performance when combining the 10 MOFA factors in a multivariate Cox regression model (assessed using cross-validation). Notably, this model yielded higher prediction accuracy than the 10 factors derived from conventional PCA or using the individual features.

Important to realize that despite the results of PCA and Factor Analysis visually look very similar, as they both provide linear variance-covariance matrix decomposition, they have quite a different math behind. While PCA is a pure matrix factorization technique which splits the total variance into orthogonal Principal Components (PCs), Factor Analysis seeks to construct hidden latent variables that generate the observed data, therefore Factor Analysis is a generative model.

 

Prepare scNMT Data Set for MOFA

In this section we will apply MOFA to the single cell multi-OMICs data set scNMT which profiles chromatine accessebility (scATACseq), DNA methylation (scBSseq) and gene expression (scRNAseq) information from the same biological cells belonging to two types: differentiating mouse embryonic stem cells (ESCs) and embrionic bodies (EBs). We will start with reading and having a look at the individual OMICs data sets:

setwd("/home/nikolay/Documents/Teaching/IntegrativeOMICs2021_ISMBECCB/workshop/session_ml/UnsupervisedOMICsIntegration")

scRNAseq<-read.delim("scRNAseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
ens2genesymbol<-read.delim("ENSEMBLE_TO_GENE_SYMBOL_MOUSE.txt")
ens2genesymbol<-ens2genesymbol[match(colnames(scRNAseq),as.character(ens2genesymbol$ensembl_gene_id)),]
colnames(scRNAseq)<-ens2genesymbol$external_gene_name
scRNAseq<-as.data.frame(t(scRNAseq))
scRNAseq[1:5,1:5]

##          ESC_A02  ESC_A03  ESC_A04  ESC_A05  ESC_A06
## Mrpl15  557.4201 130.8536 305.7356 417.1656 573.6504
## Lypla1  410.5958 138.8325 186.3237 106.0895 662.7008
## Tcea1   437.9119 356.3899 276.9121 161.8315 213.3068
## Atp6v1h 345.7199 500.5417 195.5884 133.0614 325.1376
## Oprk1     0.0000   0.0000   0.0000   0.0000   0.0000

scBSseq<-read.delim("scBSseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
scBSseq<-as.data.frame(t(scBSseq))
scBSseq[1:5,1:5]

##               ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101606896   0.00000       0       0       0       0
## 1_101935054  88.88889     100     100     100      30
## 1_101935061  88.88889     100     100     100      90
## 1_102002184 100.00000     100       0       0       0
## 1_102238204 100.00000     100     100     100     100

scATACseq<-read.delim("scATACseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
scATACseq<-as.data.frame(t(scATACseq))
scATACseq[1:5,1:5]

##             ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100392139     100     100     100     100      67
## 1_100392151     100     100     100      50      67
## 1_100668590      50      86      44      45      28
## 1_100738967      83      88      83      83      88
## 1_100994324       0       0      11      75       0

For scRNAseq OMIC layer we will only select highly expressed genes in order to remove noisy features that might contaminate the further downstream analysis. We will also perform log-transform of the data which can be seen as a mild normalization:

scRNAseq<-scRNAseq[rowMeans(scRNAseq)>=1,]
scRNAseq<-log10(scRNAseq + 1)
scRNAseq[1:5,1:5]

##          ESC_A02  ESC_A03  ESC_A04  ESC_A05  ESC_A06
## Mrpl15  2.746961 2.120092 2.486764 2.621348 2.759404
## Lypla1  2.614471 2.145608 2.272593 2.029747 2.821972
## Tcea1   2.642377 2.553142 2.443907 2.211738 2.331036
## Atp6v1h 2.539979 2.700307 2.293558 2.127304 2.513401
## Oprk1   0.000000 0.000000 0.000000 0.000000 0.000000

dim(scRNAseq)

## [1] 12145   113

Since scBSseq and scATACseq OMIC layers should be almost binary (methylated vs unmethylated for scBSseq and open vs. close for scATACseq) data sets, a good way to get rid of redundancy in the scBSseq and scATACseq data (and thus perform feature pre-selection and reduce dimensions) is to select only features with high variability:

library("mixOmics")

## Loading required package: MASS

## Loading required package: lattice

## Loading required package: ggplot2

## 
## Loaded mixOmics 6.16.3
## Thank you for using mixOmics!
## Tutorials: http://mixomics.org
## Bookdown vignette: https://mixomicsteam.github.io/Bookdown
## Questions, issues: Follow the prompts at http://mixomics.org/contact-us
## Cite us:  citation('mixOmics')

my_nearZeroVar<-nearZeroVar(as.data.frame(t(scBSseq)))
head(my_nearZeroVar$Metrics)

##             freqRatio percentUnique
## 1_101606896    112.00     1.7699115
## 1_105525627     35.00     5.3097345
## 1_105539120     26.25     5.3097345
## 1_10605572       0.00     0.8849558
## 1_107485472      0.00     0.8849558
## 1_109269874     20.00     7.9646018

dim(my_nearZeroVar$Metrics)

## [1] 3290    2

scBSseq <- scBSseq[-which(rownames(scBSseq)%in%rownames(my_nearZeroVar$Metrics)),]
scBSseq[1:5,1:5]

##               ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054  88.88889     100     100     100      30
## 1_101935061  88.88889     100     100     100      90
## 1_102002184 100.00000     100       0       0       0
## 1_102238204 100.00000     100     100     100     100
## 1_102255678 100.00000     100     100     100       0

dim(scBSseq)

## [1] 5285  113

my_nearZeroVar<-nearZeroVar(as.data.frame(t(scATACseq)), uniqueCut = 1)
head(my_nearZeroVar$Metrics)

##             freqRatio percentUnique
## 1_10401731          0     0.8849558
## 1_128314944         0     0.8849558
## 1_13628215          0     0.8849558
## 1_178804829         0     0.8849558
## 1_183332775         0     0.8849558
## 1_183944178         0     0.8849558

dim(my_nearZeroVar$Metrics)

## [1] 91  2

scATACseq <- scATACseq[-which(rownames(scATACseq)%in%rownames(my_nearZeroVar$Metrics)),]
scATACseq[1:5,1:5]

##             ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100392139     100     100     100     100      67
## 1_100392151     100     100     100      50      67
## 1_100668590      50      86      44      45      28
## 1_100738967      83      88      83      83      88
## 1_100994324       0       0      11      75       0

dim(scATACseq)

## [1] 11709   113

Let us now have a look at the histograms of individual OMICs layers in order to decide what distribution they follow and how we should model these distributions with MOFA:

hist(rowMeans(scRNAseq),breaks=100, main = "scRNAseq")

hist(rowMeans(scBSseq), breaks = 100, main = "scBSseq")

hist(rowMeans(scATACseq), breaks = 100, main = "scATACseq")

We conclude that while scRNAseq data looks fairly Gaussian (or at least exponential), we should probably model the scBSseq and scATACseq data following Bernoulli distribution as they look quite bimodal indicating the binary nature of the data, i.e. methylated vs. unmethylated for scBSseq and open vs. close for scATACseq. To make the scBSseq and scATACseq data purely Bernoulli-like, we will further make the scBSseq and scATACseq data sets binary by encoding values below 50 as 0 and above 50 as 1. Since binary data typically have vey low variation compared to continuous data, we need to remove low-variance features in this case again:

scBSseq<-ifelse(scBSseq<50,0,1)
scBSseq[1:5,1:5]

##             ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054       1       1       1       1       0
## 1_101935061       1       1       1       1       1
## 1_102002184       1       1       0       0       0
## 1_102238204       1       1       1       1       1
## 1_102255678       1       1       1       1       0

my_nearZeroVar<-nearZeroVar(as.data.frame(t(scBSseq)))
head(my_nearZeroVar$Metrics)

##             freqRatio percentUnique
## 1_110547534     21.60     1.7699115
## 1_133172191     27.25     1.7699115
## 1_139117169     27.25     1.7699115
## 1_146093078     21.60     1.7699115
## 1_194260338    112.00     1.7699115
## 1_24611678       0.00     0.8849558

dim(my_nearZeroVar$Metrics)

## [1] 761   2

scBSseq <- as.data.frame(scBSseq[-which(rownames(scBSseq)%in%rownames(my_nearZeroVar$Metrics)),])
scBSseq[1:5,1:5]

##             ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054       1       1       1       1       0
## 1_101935061       1       1       1       1       1
## 1_102002184       1       1       0       0       0
## 1_102238204       1       1       1       1       1
## 1_102255678       1       1       1       1       0

dim(scBSseq)

## [1] 4524  113

scATACseq<-ifelse(scATACseq<50,0,1)
my_nearZeroVar<-nearZeroVar(as.data.frame(t(scATACseq)))
head(my_nearZeroVar$Metrics)

##             freqRatio percentUnique
## 1_100392139    112.00      1.769912
## 1_100392151     27.25      1.769912
## 1_100738967     55.50      1.769912
## 1_101897864     55.50      1.769912
## 1_102283335     55.50      1.769912
## 1_102563492    112.00      1.769912

dim(my_nearZeroVar$Metrics)

## [1] 2238    2

scATACseq <- as.data.frame(scATACseq[-which(rownames(scATACseq)%in%rownames(my_nearZeroVar$Metrics)),])
scATACseq[1:5,1:5]

##             ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100668590       1       1       0       0       0
## 1_100994324       0       0       0       1       0
## 1_100994328       0       0       0       0       0
## 1_100994336       0       0       0       0       0
## 1_100994343       0       0       0       0       1

dim(scATACseq)

## [1] 9471  113

Now data cleaning step is finished and the OMICs are prepared to be integrated in unsupervised way with MOFA.

Run MOFA on scNMT Data Set

Let us continue with creating MOFA object. MOFA allows for a handy overview of the data by displaying dimensions of each OMIC.

library("MOFA2")

## 
## Attaching package: 'MOFA2'

## The following object is masked from 'package:stats':
## 
##     predict

omics<-list(scRNAseq=as.matrix(scRNAseq),scBSseq=as.matrix(scBSseq),scATACseq=as.matrix(scATACseq))
lapply(omics,dim)

## $scRNAseq
## [1] 12145   113
## 
## $scBSseq
## [1] 4524  113
## 
## $scATACseq
## [1] 9471  113

MOFAobject <- create_mofa_from_matrix(omics)
plot_data_overview(MOFAobject, )

Now we will use some model training parameters such as distribution types corresponding to each OMIC, number of iterations, number of factors to be computed etc.

#DEFINE DATA OPTIONS
DataOptions <- get_default_data_options(MOFAobject)
DataOptions

## $scale_views
## [1] FALSE
## 
## $scale_groups
## [1] FALSE
## 
## $center_groups
## [1] TRUE
## 
## $use_float32
## [1] FALSE
## 
## $views
## [1] "scRNAseq"  "scBSseq"   "scATACseq"
## 
## $groups
## [1] "group1"

#DEFINE MODEL OPTIONS
ModelOptions <- get_default_model_options(MOFAobject)
mydistr<-c("gaussian","bernoulli","bernoulli")
names(mydistr)<-c("scRNAseq","scBSseq","scATACseq")
ModelOptions$likelihoods<-mydistr
ModelOptions$num_factors<-20
ModelOptions

## $likelihoods
##    scRNAseq     scBSseq   scATACseq 
##  "gaussian" "bernoulli" "bernoulli" 
## 
## $num_factors
## [1] 20
## 
## $spikeslab_factors
## [1] FALSE
## 
## $spikeslab_weights
## [1] TRUE
## 
## $ard_factors
## [1] FALSE
## 
## $ard_weights
## [1] TRUE

#DEFINE TRAIN OPTIONS
TrainOptions <- get_default_training_options(MOFAobject)
TrainOptions$seed <- 2018
# Automatically drop factors that explain less than 3% of variance in all omics
TrainOptions$drop_factor_threshold <- 0.03
# TrainOptions$tolerance <- 0.1
TrainOptions$maxiter<-1000
TrainOptions$verbose<-TRUE
TrainOptions

## $maxiter
## [1] 1000
## 
## $convergence_mode
## [1] "fast"
## 
## $drop_factor_threshold
## [1] 0.03
## 
## $verbose
## [1] TRUE
## 
## $startELBO
## [1] 1
## 
## $freqELBO
## [1] 5
## 
## $stochastic
## [1] FALSE
## 
## $gpu_mode
## [1] FALSE
## 
## $seed
## [1] 2018
## 
## $outfile
## NULL
## 
## $weight_views
## [1] FALSE
## 
## $save_interrupted
## [1] FALSE

Finally, we are ready to run MOFA:

MOFAobject <- prepare_mofa(MOFAobject, data_options = DataOptions, model_options = ModelOptions, training_options = TrainOptions)

## Checking data options...

## Checking training options...

## Checking model options...

#reticulate::use_python("/home/nikolay/miniconda3/envs/HOPS/bin/python3")
#reticulate::use_condaenv("/home/nikolay/miniconda3/envs/HOPS", required = TRUE)
MOFAobject <- run_mofa(MOFAobject, use_basilisk = FALSE)

## Connecting to the mofapy2 python package using reticulate (use_basilisk = FALSE)... 
##     Please make sure to manually specify the right python binary when loading R with reticulate::use_python(..., force=TRUE) or the right conda environment with reticulate::use_condaenv(..., force=TRUE)
##     If you prefer to let us automatically install a conda environment with 'mofapy2' installed using the 'basilisk' package, please use the argument 'use_basilisk = TRUE'

MOFAobject

## Trained MOFA with the following characteristics: 
##  Number of views: 3 
##  Views names: scRNAseq scBSseq scATACseq 
##  Number of features (per view): 12145 4524 9471 
##  Number of groups: 1 
##  Groups names: group1 
##  Number of samples (per group): 113 
##  Number of factors: 3

Let us check the covariance of the OMICs layers:

#ANALYZE RESULTS OF MOFA INTEGRATION
head(get_variance_explained(MOFAobject)$r2_total[[1]])

##  scRNAseq   scBSseq scATACseq 
## 13.297925  5.474570  4.624678

head(get_variance_explained(MOFAobject)$r2_per_factor[[1]])

##         scRNAseq    scBSseq  scATACseq
## Factor1 7.008601 5.41874498 4.46349678
## Factor2 3.547561 0.02761408 0.12667731
## Factor3 3.163781 0.02591241 0.03993401

MOFAobject@cache$variance_explained

## $r2_total
## $r2_total$group1
##  scRNAseq   scBSseq scATACseq 
## 13.297925  5.474570  4.624678 
## 
## 
## $r2_per_factor
## $r2_per_factor$group1
##         scRNAseq    scBSseq  scATACseq
## Factor1 7.008601 5.41874498 4.46349678
## Factor2 3.547561 0.02761408 0.12667731
## Factor3 3.163781 0.02591241 0.03993401

plot_list=plot_variance_explained(MOFAobject, x='view', y='factor', plot_total = T)
plot_list[[2]] #total variance

plot_list[[1]] #variance by factor

We can see that scRNAseq provides the largest variation (13%) in the total integrative OMICs scNMT data set, scBSseq and scATACseq contribute around 5% of variation. We also observe that MOFA selected 3 Latent Factors (LFs), scRNAseq contributes to all of them while scBSseq and scATACseq contributes only to the first LF. Now let us interpret the LFs by displaying the features associated with each LF:

NumFactors<-dim(get_factors(MOFAobject)$group1)[2]

plot_weights(MOFAobject, view = "scRNAseq", factor = 1)

plot_top_weights(object = MOFAobject, view = "scRNAseq", factor = 1, nfeatures = 10)

plot_data_heatmap(object = MOFAobject, view = "scRNAseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)

plot_weights(MOFAobject, view = "scBSseq", factor = 1)

plot_top_weights(object = MOFAobject, view = "scBSseq", factor = 1, nfeatures = 10)

plot_data_heatmap(object = MOFAobject, view = "scBSseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)

plot_weights(MOFAobject, view = "scATACseq", factor = 1)

plot_top_weights(object = MOFAobject, view = "scATACseq", factor = 1, nfeatures = 10)

plot_data_heatmap(object = MOFAobject, view = "scATACseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)

Above we ranked by their importance and visualized features from each OMIC that contribute the most to the Factor 1 where all three OMICs (scRNAseq, scBSseq and scATACseq) are strongly correlated. Therefore one can expect those features to have something in common, i.e. they should also be correlated or involved in the same biological pathway.

Let us also display the cells in the low-dimensional latent space. for this purpose we will input MOFA factors into tSNE and UMAP and display the cells as a consensus plot for scRNAseq, scBSseq and scATACseq:

controls<-c("EB_P1D12","EB_P1E12","EB_P1F12","EB_P1G12","EB_P2B12","EB_P2D12","EB_P2E12","EB_P2F12","EB_P2G12",
            "ESC_B12","ESC_C12","ESC_D12")
colnames(scRNAseq)<-ifelse(colnames(scRNAseq)%in%controls,"CONTROL_CELL",colnames(scRNAseq))
cell_types<-matrix(unlist(strsplit(colnames(scRNAseq),"_")),ncol=2,byrow=TRUE)[,1]

sample_metadata <- data.frame(sample = samples_names(MOFAobject)[[1]],condition = cell_types)
samples_metadata(MOFAobject) <- sample_metadata
head(samples_metadata(MOFAobject), n=3)

##    sample condition  group
## 1 ESC_A02       ESC group1
## 2 ESC_A03       ESC group1
## 3 ESC_A04       ESC group1

set.seed(12345)
MOFAobject <- run_umap(MOFAobject, n_neighbors = round(sqrt(length(cell_types)),0), min_dist = 1)
MOFAobject <- run_tsne(MOFAobject, perplexity = round(sqrt(length(cell_types)),0))

plot_dimred(MOFAobject, method = "TSNE", dot_size = 3, color_by = "condition")

plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "condition")

#plot(MOFAobject@dim_red$TSNE$TSNE2~MOFAobject@dim_red$TSNE$TSNE1, col=as.factor(cell_types), xlab = "TSNE1", ylab = "TSNE2", pch=19,
#     main="tSNE representation of MOFA integration of scRNAseq, scBSseq and scATACseq", cex = 2)
#legend("bottomleft", inset=.02, levels(as.factor(cell_types)), fill=c("black","red","green"), cex=1.5)
#plot(MOFAobject@dim_red$UMAP$UMAP2~MOFAobject@dim_red$UMAP$UMAP1, col=as.factor(cell_types), xlab = "UMAP1", ylab = "UMAP2", pch=19,
#     main="UMAP representation of MOFA integration of scRNAseq, scBSseq and scATACseq", cex = 2)
#legend("bottomleft", inset=.02, levels(as.factor(cell_types)), fill=c("black","red","green"), cex=1.5)

We conclude that the MOFA unsupervised integrative OMICs approach can distinguish between all three types of cells: ESCs, EBs and Controls. Now let us see what MOFA factors contribute to different clusters:

plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor1")

plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor2")

plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor3")

We can see that the ESC cluster is almost entirely driven by the MOFA Factor 1 meaning that all three OMICs (scRNAseq, scBSseq and scATACseq) contribute more or less equally into this cluster. In contrast the Control cluster is driven by the factor 2, which is mainly due to scRNAseq contribution. Finally, Factor 3 contributes more or less equally to all the three clusters implying that scRNAseq impacts all the three clusters.

Session Info

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=sv_SE.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=sv_SE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=sv_SE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=sv_SE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] MOFA2_1.3.4     mixOmics_6.16.3 ggplot2_3.3.5   lattice_0.20-44
## [5] MASS_7.3-54    
## 
## loaded via a namespace (and not attached):
##  [1] MatrixGenerics_1.4.3 sass_0.4.0           tidyr_1.1.3         
##  [4] jsonlite_1.7.2       ellipse_0.4.2        bslib_0.3.0         
##  [7] assertthat_0.2.1     highr_0.9            stats4_4.1.0        
## [10] yaml_2.2.1           ggrepel_0.9.1        corrplot_0.90       
## [13] pillar_1.6.2         glue_1.4.2           reticulate_1.20     
## [16] digest_0.6.27        RColorBrewer_1.1-2   colorspace_2.0-2    
## [19] cowplot_1.1.1        htmltools_0.5.2      Matrix_1.3-4        
## [22] plyr_1.8.6           pkgconfig_2.0.3      pheatmap_1.0.12     
## [25] dir.expiry_1.0.0     purrr_0.3.4          corpcor_1.6.9       
## [28] scales_1.1.1         HDF5Array_1.20.0     RSpectra_0.16-0     
## [31] Rtsne_0.15           BiocParallel_1.26.2  tibble_3.1.4        
## [34] farver_2.1.0         generics_0.1.0       IRanges_2.26.0      
## [37] ellipsis_0.3.2       withr_2.4.2          BiocGenerics_0.38.0 
## [40] magrittr_2.0.1       crayon_1.4.1         evaluate_0.14       
## [43] fansi_0.5.0          forcats_0.5.1        tools_4.1.0         
## [46] lifecycle_1.0.0      matrixStats_0.60.1   basilisk.utils_1.4.0
## [49] stringr_1.4.0        Rhdf5lib_1.14.2      S4Vectors_0.30.0    
## [52] munsell_0.5.0        DelayedArray_0.18.0  compiler_4.1.0      
## [55] jquerylib_0.1.4      rlang_0.4.11         rhdf5_2.36.0        
## [58] grid_4.1.0           rhdf5filters_1.4.0   RcppAnnoy_0.0.19    
## [61] igraph_1.2.6         labeling_0.4.2       rmarkdown_2.10      
## [64] basilisk_1.4.0       codetools_0.2-18     gtable_0.3.0        
## [67] DBI_1.1.1            rARPACK_0.11-0       reshape2_1.4.4      
## [70] R6_2.5.1             gridExtra_2.3        knitr_1.33          
## [73] dplyr_1.0.7          uwot_0.1.10          fastmap_1.1.0       
## [76] utf8_1.2.2           filelock_1.0.2       stringi_1.7.4       
## [79] parallel_4.1.0       Rcpp_1.0.7           vctrs_0.3.8         
## [82] png_0.1-7            tidyselect_1.1.1     xfun_0.25